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ABSTRACT 

We model the evolution of infrared galaxies using a phenomenological approach to 
match the observed source counts at different infrared wavelengths. In order to do 
that, we introduce a new algorithm for reproducing source counts which is based on 
direct integration of probability distributions rather than using Monte-Carlo sampling. 
We also construct a simple model for the evolution of the luminosity function and the 
colour distribution of infrared galaxies which utilizes a minimum number of free pa- 
rameters; we analyze how each of these parameters is constrained by observational 
data. The model is based on pure luminosity evolution, adopts the Dale & Helou 
Spectral Energy Distribution (SED) templates, allowing for evolution in the infrared 
luminosity-colour relation, but does not consider Active Galactic Nuclei as a separately 
evolving population. We find that the 850 /im source counts and redshift distribution 
depend strongly on the shape of the luminosity evolution function, but only weakly on 
the details of the SEDs. Based on this observation, we derive the best-fit evolutionary 
model using the 850 /zm counts and redshift distribution as constraints. Because of the 
strong negative ii'-correction combined with the strong luminosity evolution, the fit 
has considerable sensitivity to the sub-L* population at high redshift, and our best-fit 
shows a flattening of the faint end of the luminosity function towards high redshifts. 
Furthermore, our best-fit model requires a colour evolution which implies the typical 
dust temperatures of objects with the same luminosities to decrease with redshift. We 
then compare our best-fit model to observed source counts at shorter and longer wave- 
lengths which indicates our model reproduces the 70 /im and 1100 /xm source counts 
remarkably well, but under-produces the counts at intermediate wavelengths. Analy- 
sis reveals that the discrepancy arises at low redshifts, indicating that revision of the 
adopted SED library towards lower dust temperatures (at a fixed infrared luminosity) 
is required. This modification is equivalent to a population of cold galaxies existing at 
low redshifts, as also indicated by recent Herschel results, which are underrepresented 
in IRAS sample. We show that the modified model successfully reproduces the source 
counts in a wide range of IR and submm wavelengths. 

Key words: galaxies: evolution - galaxies: general - galaxies: starburst - infrared: 
galaxies - submillimetre: galaxies 



1 INTRODUCTION 

Although dust is an unimportant component in the mass 
budget of galaxies, its presence radically alters the emer- 
gent spectrum of star forming galaxies. Since stars are born 
in dusty clouds, most of the energy of young stars is absorbed 
by dust particles, which are heated by the absorption pro- 
cess and radiate their energy away by thermal emission at 
infrared (IR) and submillimetre (submm) wavelengths. As 
a result, the infrared radiation of star forming galaxies is a 
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useful measure of the massive star formation rate. Mainly 
due to the atmospheric opacity, the thermal radiation from 
dust could not be systematically studied until the launch 
of IRAS in the mid-1980s, which provided the first compre- 
hensive view of the nearby dusty Universe. Ten years later 
COBE discovered the Cosmic Infrar ed Background (CIB) 
(|Puget et al.ll 19961 : iFixsen et aniigGSl ) and it turned out the 
observed power of the CIB is comparable to what can be de- 
duced from the optical Universe. This was in contrast to the 
observations of local galaxies which suggest o nly one third of 
the e nergy output of galaxies is in IR bands IjLagache et al.l 
[2OO5I). Moreover, the relatively fiat slope of the CIB at long 



2 A. Rahmati and P.P. van der Werf 



wavelengths indicated a population of dusty galaxies which 
are d istributed over a wide range of redshifts IjGispert et alj 
I2OO0I). Thanks to various large surveys performed with dif- 
ferent satellites and ground based observatories, this back- 
ground radiation is now partly resolved into point sources 
at different wavelengths. While at shorter wavelengths the 
emission is mainly coming from local and low redshift galax- 
ies, longer wavelengths contain information about larger dis- 
tances. 

The shape of emission spectrum of warm dust parti- 
cles resembles a modified blackbody spectrum with a peak 
varying with the typical dust temperature which is ob- 
served to be around T ~ 30 — 40K in galaxies; therefore 
the far-IR (FIR) and submm spectrum of a typical dusty 
galaxy consists of a peak at rest-frame wavelengths around 
A ~ 100 — 200/^m which drops on both sides (see Figure [T]). 
While the presence of other emitters like polycyclic aromatic 
hydrocarbons (PAHS) and Active Galactic Nuclei (AGNs) 
complicates the shape of spectra at shorter wavelengths, at 
submm wavelengths the spectra behave like modified black- 
bodies and their amplitudes drop steeply. In fact, because 
of this steep falloff in the Rayleigh- Jeans tail of the Spectral 
Energy Distribution (SED), the observed flux density at a 
fixed submm wavelength can rise by moving the SED in red- 
shift space (the so called AT-correction) , which counteracts 
the cosmological dimming. It turns out that the interplay 
between these two processes enables us to observe galaxies 
at submm wavelengths out to very high redshifts (see also 
the discussion in Section [S] and Figure [7|. 

After the first observations of SCUBA at 850 /^m con- 
firmed the import ance of submm galaxies at high redshifts 
l|SmaiI et al.ll 19971 ) , there were many subsequent surveys us- 
ing differe nt instruments wh i ch explored different cosmolog- 
ical fields (Smail ot al."2002'; 'Webla ct al."2003'; 'Borvs et al.l 
2003; Grcvc ct al. 2004; Coppin ct al. 2006; Bcrtoldi ct~aL| 



20071: lAustermann et aklboiol : IScott et al.ll201Gl : IVieira et al 



2010|) and also some surveys which used gravitational lens- 
ing techniques to extend the observable submm Universe 
to sub-mjy fl uxes (Small ct al. 1997; Cha pman et al.ll2002l : 
iKnudsen et a l. 2008; Johansson ct al. 201l|). 

While low angular resolution makes individual iden- 
tifications and spectroscopy a daunting task, the surface 
density of sources as a function of brightness (i.e., the 
source counts) can be readily analysed and contains signifi- 
cant information about the population properties and their 
evolution. One can assume simple smooth parametrized 
models for the evolution of dusty galaxies and relate 
their low redshift observed properties (e.g., the total IR 
luminosity function, which we simply call the luminos- 
ity func tion, or in short LF h e reafter) to their source 
counts (JBlain fc LongaiJ Il993l: JGuiderdoni et al.l [l997l; 
Blain et all |l 999'; 'Chary fc Elb aj I2OOII : iRowan-Robinson' 



2001 



2005 



Dole et al. 2003; Lagache et al. 2004; Lewis ct al 



Le Borgne et all l2009l : IValiante et all |2009| ). Such 



backward evolution models usually combine SED templates 
and the low redshift properties of IR galaxies which are al- 
lowed to change with redshift based on an assumed para- 
metric form, together with Monte-Carlo techniques to re- 
produce the observed source counts. These models are 
therefore purely phenomenological, and suppress the un- 
derlying physics. Their power lies exclusively in providing 
a parametrized description of statistical properties of the 



galaxy population under study, and the evolution of these 
properties. The results of such modeling provides a descrip- 
tion of the constraints that must be satisfied (in a statistical 
sense) and could be explored further by more physically mo- 
tivated models such as hydrodynamical simulations embed- 
ded in a ACDM cosmology, with subgrid presc r iption s for 
star formation and feedback (e.g. ISchave et al.l l|20ld )). It 
is important to realize that backward evolution models are 
constrained only over limited observational parameter space 
(e.g., at particular observing wavelengths, flux levels, red- 
shift intervals, etc.), and are ignorant of the laws of physics. 
Therefore it is dangerous to use them outside of their es- 
tablished validity ranges; in other words, these models have 
descriptive power but little explanatory power. As a result, 
one of the most instructive uses of these models is in analyz- 
ing where they fail to reproduce the observations, and how 
they can be modified to correct for these failures which im- 
plies that the underlying assumptions must be revised. This 
is the approach that we use in the present paper. 

Since nowadays many and various observational con- 
straints exist (e.g. source counts at various IR and submm 
wavelengths, colour distributions, redshift distribution), 
backward evolution models can reach considerable levels 
of sophistication. However, early backward evolution mod- 
els used only a few SED templates (or sometimes even 
only one SED) to represent the whole population of dusty 
galaxies. This approach neglects the fact that dusty galax- 
ies are not identical and span a variety of dust tempera- 
tures and hence SED shapes. As a result, such models can- 
not probe the possible evolution in the SED properties of 
submm galaxies, for which increasing obser vational support 



has been found d uring the l ast few years (I Chapman et al 



20051 : iPope et al.|[200 6: Chapi n et al.ll2009l:ISvmeonidis et al 
2OO9I : ISevmour et al.ii201Q : iHwang et al.ll2010l ') 



In addition, existing backward evolution models typ- 
ically use only "luminosity" and/or "density" evolution 
which respectively evolves the characteristic luminosity (i.e., 
I/*) and the amplitude of the LF without changing its shape. 
In other words, they do not consider any evolution in the 
shape of the luminosity function. Finally, the intrinsically 
slow Monte-Carlo methods used in many of the models make 
the search of parameter space for the best-fit model labori- 
ous and tedious. 

In this paper, we use a new and fast algorithm which is 
different from conventional Monte-Carlo based algorithms, 
for calculating the source counts for a given backward evo- 
lution model. We also use a complete library of IR SED 
templates which form a representative sample of observed 
galaxies, at least at low redshifts. This will enable us to ad- 
dress questions about the evolution of colour distribution of 
dusty galaxies during the history of the Universe. We also 
consider a new evolution form for the LF which allows us 
to constrain evolution in the shape of the LF in addition to 
quantifying its amplitude changes. 

The structure of the paper is as follows: in fj2]we present 
different ingredients of our parametric colour-luminosity 
function (CLF) evolution model and introduce a new algo- 
rithm for the fast calculation of the source counts for a given 
CLF. Then, after choosing our observational constraints for 
850 ^m objects in ^we try to find a model consistent with 
those constraints by studying the sensitivity of the produced 
source counts to different parameters in SQ] After finding a 
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850 ^m-constrained best-fit model, we test its performance 
in producing source counts at otiier wavelengths in ^which 
leads us to a model capable of producing source counts in a 
wide range of IR and submm wavelengths. Finally, after dis- 
cussing the astrophysical implications of our best-fit model 
in 3ni we end the paper with concluding remarks. Through- 
out this paper, we use the standard ACDM cosmology with 
the parameters ilm = 0.3, Qa = 0.7 and h = 0.75. 



2 MODEL INGREDIENTS 

Although luminosity is the main parameter to differentiate 
between different galaxies, it is an integrated property of 
SED and cannot distinguish between different SED shapes 
corresponding to different physical processes taking place 
in different objects with similar outgoing integrated energy 
fluxes. U sing colour indicators can resolve this degeneracy. 
iDale et a l. (2001) demonstrated that the i?(60, 100) which 
is defined as log(S60/jm/S'iooMm)! is the best single param- 
eter characterization of IRAS galaxies. Moreover, it was 
shown that IRAS galaxies exhibit a slowly varying correla- 
tion between R{60, 100) and luminosity such that objects 
with larger charac teristic luminosit i es have warmer char- 
acteristic colours (iDale et al.l l200ll : IChapman et al.l l2005l : 
IChapin et al.l 12009 ). Based on these facts, we construct a 
model for the Colour-Luminosity Function (hereafter CLE). 
This model consists of the observed CLE in the local Uni- 
verse and a parametric evolution function. Then, by adopt- 
ing a suitable set of SED models, we calculate the source 
count of IR objects at different wavelengths and compare 
the results with observations. 



2.1 CLF at z = 

The local volume density of IRAS galaxies, ^o{L,C), can 
be parametrized as a function of total IR luminosity, L, and 
R{60, 100) colour, C, where the total infrared luminosity is 
calculated by int egrating over the SED from 3 to 1100 /im 
([Dale et al.ll200ll ). Furthermore, it is possible to express the 
local colour-luminosity distribution as the product of the 
local luminosity function, $o(i), and the local conditional 
probability of a galaxy having the colour C given the lumi- 
nosity L, Po{C\L), 



ML,C) = <!>o{L)Po{C\L). 



(1) 



Since the IRAS galaxies represent an almost com- 
plete sam£le__of_IR__galaxies out to the redshift ~ 
0.1 ([Saunders et al.l Il990h , we use the luminosity and 
colour functions found by analyzing a ff ux limited sam- 
ple o f Seoum > 1 . 2Jy I R AS galaxies (^Chap man et al.l 
20031: IChapin et al.1 |2009| ). [Chapman et al. (2001) and 



Chapin et al.l ( 2003 ) analyze this sample which covers most 
of the sky, using an accessible volume technique for finding 
the LF and fit a dual power law function to the observed 
luminosity distributio n. The parametr ic form of luminosity 
function based oU iChapin et al.l (|2009l ) is given by 



$o(i) = P*( 



L,' 



"(1+1^)-^ 



(2) 



characterize the power-laws at the faint {L < Lt) and bright 
[L > L«) ends, respe c tively . 

IChapman et al.1 (|2003l ) and IChapin et ahl (12009') also 
found a Gaussian representation for the colour distribution 
of IRAS galaxies 

Po(C7!L) = -=^exphix(^^^)^ (3) 

V27rcrc ^ ctc 



where L, = 5.14 x W^'^Lq is the characteristic knee luminos- 
ity, p* = 1.22 X 10~^^Mpc~"^Lq^ is the number density nor- 
malization of the function at L,, and a = 2.59 and /3 = 2.65 



where Co can be represented by a dual power law function 

Co = a-51og(l + :^) + 7log(l-f -^), (4) 

withC. =-0.48, (5 = -0.06,7 = 0.21 and L' = 3.2x10^Lq, 
and the distribution width, CTc, is expressed as 

ac = af(l - 2-^'/^) + ab(l - 2"^/^'), (5) 

where erf = 0.2 and at = 0.128 (|Chapin et al.ll2009l ). 



2.2 CLF evolution 

The CLF introduced in the previous section, contains in- 
formation about the distribution of IR galaxies only in the 
nearby Universe. For modeling the IR Universe at higher 
redshifts, its evolution must be modeled. The necessity of 
the CLF evolution with redshift is shown directly by the 
fact that the observed power of the the CIB is comparable 
to what can be deduced from the optical cosmic background 
cannot be explained by our understanding of the local Uni- 
verse which indicates the infrared output of galaxies is only 
one third of their optical output (Laga che et al.ll2005i ). Fur- 
thermore, several studies have shown for a fixed total IR lu- 
minosity the typical temperature of infrared sources is lower 
at higher r edshifts (Chapman et al. 2005; Pope et al .1120061: 
Chapin et al. 2009; Symeonidis et al. 2009; S eymo ur et al.l 
20101 : iHwang et al.ll2O10l : lAmblard et al.ll2010l ). which advo- 
cates an IR colour evolution. Therefore, we need to adopt 
a reasonable form of evolution in the luminosity and colour 
distributions to reproduce correctly the CLF evolution of IR 
galaxies with redshift. 

There are three different ways to evolve the local lu- 
minosity function of IR galaxies: (i) changing p« with red- 
shift (i.e., density evolution); (ii) changing L* with redshift 
(i.e., luminosity evolution) and (iii) changing the bright and 
faint end slopes (i.e., a and /3) with redshift. While the 
density and luminosity evolution change the abundances of 
all sources independent of their luminosities, they leave the 
shape of the LF unchanged. Such models are called transla- 
tional models since they amount to only a translation of the 
LF in parameter space. Going beyond translational models, 
a variation of the shape of LF with redshift (i.e., varying 
a and /3 slopes in equation [Sjl can be used to change the 
relative contributions of bright and faint sources at different 
redshifts. 

As discussed bv lBlain et al.l (|l999r ). luminosity and den- 
sity evolution affect the CIB predicted by the models simi- 
larly, but luminosity evolution has a much stronger effect on 
the source counts. Combining source counts and integrated 
background can therefore distinguish between luminosity 
and density evolution. The result of this analysis shows 
luminosity evolution must strongly dominate over density 
evolution, since pure density evolution consistent with the 
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observed 850 ^m source counts would ov erpredict the inte - 
grated background by a factor 50 to 100 IjBlain et al.ll 19991 ). 
We therefore assume negligible density evolution. However, 
as we will show in Section 4, luminosity evolution is not suf- 
ficient to reproduce the correct source count and redshift 
distribution of submm galaxies, which forces us to drop the 
assumption of a purely translational LF evolution model. 
Moreover, there is no reason to assume that the slopes of 
LF at faint and bright ends remain the same at all redshifts. 
Therefore, we allow them to change in our model and intro- 
duce a redshift dependent LF which can be written as 

.^iL,z) = p4^^)'--^{l + ^^)-^% (6) 

g{z)L* g{z)L, 

which is similar to equation ((2|) but now q^ and [3^ are chang- 
ing with redshift and L* is multiplied by g{z), which we refer 
to as the luminosity evolution function. 

Since the average total-IR-luminosity density and its as- 
sociated star formation density closely follow the luminosity 
evolution function, we choose a form of luminosity evolution 
which is similar to the observed evolution o f the cosmic star 
formation history (see Hopkins & Beacon (120061 )): in other 
words, we assume that the luminosity evolution is a func- 
tion which is increasing at low redshifts and after reaching 
its maximum turns into a decreasing function of look-back 
time at high redshifts: 



9{z) = 




Zb 



if Z < Za 

if Za < Z < Zt 

if 26 < Z 



(7) 



where n, m, Za and zj, are constants. As we will show in 
Section 14.21 the source count is not strongly sensitive to 
the model properties at high redshifts. In other words, the 
sensitivity of the source count calculation to the difference 
between Za and zt is much less than its high sensitivity to 
the value of Za itself. Moreover, the model outputs are not 
highly sensitive to the slope of g{z) at high redshifts if it 
remains negative (i.e., m < 0). Therefore, we can simplify 
the model by assigning suitable fixed values to Zb — Za and m, 
without any significant change in its fiexibility. We choose 
Zb — Za = ^ and m = —1, knowing that any different choice 
for those values can be compensated by a very small change 
in Za or/and n. we will discus this issue in more detail in 
Section 1321 

We also adopt linear forms for changing a^ and (3^ with 
redshift 



Oz — 2.59 + aaZ, 



Pz = 2.65 -I- a/jz, 



(8) 
(9) 



where a^ and a/s are constants. We apply the slope evolution 
only for z < '"2"'' noting that the observed source count 
is not very sensitive to the exact properties of our model at 
very high redshifts and extending the evolution of LF slopes 
to even higher redshifts does not change our results. 

As mentioned earlier, some studies have shown that 
high redshift IR galaxies have lower dust temperatures than 
low redshift galaxies, at a fixed total IR luminosity. In our 
model, in order to evolve the c olour distributi o n acc ordingly, 
we adopt a relation similar to lValiante et al.l l)2009l l to shift 
the centre of the colour distribution for a given luminos- 
ity towards colours associated with lower dust temperatures 



03 




10 100 

X(/Lim) 



1000 



Figure 1. Some S EP samples from the template library of 
IPale fc Helod l|2002h which we use in our model. The blue, green, 
orange and red lines (from top to bottom) are respectively asso- 
ciated with colours C = 0.16, 0.04, -0.27 and -0.54. The nor- 
malization is arbitrary. 



(i.e., smaller _R(60, 100) values) at higher redshifts (see Fig- 
ure lllfl . In this way, the colour distribution of IR galaxies 
can be written as 



P{C I L) 



1 . 1 ,C- 

„^exp[-- X (— - 
27r(Tc 2 a, 



Co. 



where a^ is given by equation ([5]) and Cq is 

Co = a - 5iog(i + ^ I ' ) + 7iog(i + 



L 



L'{l + z) 



(10) 



:l 



(11) 

where ui is a constant. Alth ough for colour evolu tion we 
adopt a similar formalism to IValiante et al.l ((20091), unlike 
them we allow w to vary as a free parameter which enables 
us to constrain the colour evolution as well. 

Based on the above formalism, the general CLF of IR 
galaxies can be written as 

<i>{L,C,z) = <i>{L,z)P{C\L), (12) 

where ^{L, z) and P{C\L) are defined by equations (|5))- pip . 

2.3 SED Model 

Several studies have shown that AGNs do not domin ate the 
FIR e nergy output o f the U niverse ( Swinban k et al.l 



2004 [Alexander et ai] l2005l: IChapman et al.l 
Lutz et al.l l2005l: IValiante et all l2007l: iPope et al.l 




Menendez-Delmestre et al.l l2009l : iFadda et all I201C 
J auzac et al.l 1201(1 ). Therefore, in order to avoid un- 
necessary complexity in our model we adopt a single 
population of galaxies, modeled with only one family of 
SEDs representing star forming galaxies. 

This choice is justified for long wavelengths, for instance 
at 850 lira, by considering the fact that at those wavelengths 
the AGN contribution to the observed fiux is negligible. In 
other words, to be able to change the 850 ^m fiuxes and 
source counts, AGN should be the dominant contributor to 
the SED at rest-frame wavelengths longer than ~ 200 pm 
which is highly unlikely. 
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On the other hand, excluding AGN contribution could 
be potentially important at shorter wavelengths: if we adopt 
a simple assumption where the AGN cont inuum is weU rep- 
resented by a simp le torus model (Efs tathiou et al.lll995l : 
IValiante et al.ll2009r ). then a strong enough AGN could cre- 
ate a bump, on top of the starburst SED, at rest-frame wave- 
length ranges ~ 10 — 5 lira (see Figure 2 in Efstathiou et alj 
l|l995[l and Figure 9 in lValiante et all (|2009 l')). This feature, 
together with increasing AGN luminosity with redshift could 
modify the observed source counts at observed wavelengths 
shorter than 200 iim but still is unlikely to affect submm 
counts. 

However. iMulIanev et al.l (|201lh recently used the deep- 
est available Herschel survey and showed that for a sample 
of X-ray selected AGNs up to z ~ 3, the observed 100 /im 
and 160 /^m fluxes are not contaminated by AGN. Conse- 
quently, even if all the galaxies which contribute to the FIR 
and submm source counts host AGN, their observed fluxes 
is driven by star formation activity at wavelengths around 
100 nm or longer. 

In order to get ac curate SED templates for star-forming 
lalaxi es, we use the iDale et al.l (|200ll ) and [Dale & Hclou 
2002f) SED models which are produced by a semi-empirical 
method to represent spectral energy distributi on of star- 
formin g galaxies in the IR region of spectrum. I Dale et al.l 
(2001) add up emission profiles of different dust families 
(i.e. large grains, very small grains and polycyclic aromatic 
hydrocarbons) which are exposed to a range of radiation 
field strengths in a parametrized form to generate the SED 
of different star-forming systems. In this way, it is possi- 
ble to produce the SED corresponding to each _R(60, 100) 
colour and scale it to the desired total infrared luminos- 
ity. In our model, we use spectral templates taken from the 
I Dale fc Heloul ([2003) catalog which provides 64 normalized 
SEDs with different _R(60, 100) colours, ranging from —0.54 
to 0.21. Some SED examples taken from this template set 
are illustrated in Figure [l] 



frame luminosity density (WHz~^) of the object with total 
luminosity L and colour C at wavelength Ao ~ Xohs/{z + 1). 
Now consider a cell defined by redshift interval Az — 
Z2 — z\ (where z\ < 22), luminosity interval AL — L2 — Li 
and colour interval AC — C2 — Ci. Assuming negligible 
colour and luminosity evolution between zi and 22, the av- 
erage detection probability in the cell is equal to the fraction 
of detectable objects in that cell which is: 



1 

D^(z„ 



.)-D-'(^l) 



D3(z2)-D^{zi) 







if Zmax > Z2 

if Zl < 2max < Z2 

n -2'max ^ Zl 



(15) 



where D{z) is the proper distance equivalent to redshift z. 
For writing equation 1151 we assumed galaxies to be dis- 
tributed uniformly in space between redshifts 21 and 22. 
Consequently, based on equation 1131 the number of objects 
which are contributing to the source count in each cell of 
the redshift-colour-luminosity space is 



AN{> Stb) = Q X $(L,C,2)ALACA1/, 



(16) 



where AV is the volume corresponding to the redshift inter- 
val A2 and $(!/, C, 2) is the average CLE in the cell, which 
for small enough AL, AC and A2 can be written as 



$(L,C,2) = <1.(L,C,2) 



(17) 



Then, the total source count is obtained by summing over 
the contribution from all cells (see Appendix A). 

As a result of using Q and its special form as expressed 
in equation ll5l our algorithm computes the source count for 
continuously distributed sources in the Universe, indepen- 
dent of the size of A2 which is used in equation [TS] Thanks 
to this property, the computational cost required for the 
source count calculation is reduced significantly, and the 
small size of A2 is only necessary for accurate calculation 
of ii'-corrected SEDs and CLFs at different redshifts and 
not to guarantee a uniform distribution of sources (see also 
Appendix A). 



2.4 The algorithm 

Given the distribution of objects in the Universe, and their 
associated SEDs, it is possible to calculate the number of 
sources which have observed fluxes above some detection 
threshold, Sth, at a given wavelength A — Aobs: 



N{> Stir) = 



Q X $(L, C, z)—dzdLdC, (13) 
az 



where V is the volume and Q is the probability that a 
source with luminosity L, colour C and redshift 2 has an 
observed flux density greater than Sth at that wavelength 
(i.e., its detection probability). At each point of redshift- 
colour-luminosity space, Q is either 1 or for any particular 
galaxy, but it is useful to think of Q as the average proba- 
bility of detection for galaxies in each cell of that space. 

In order to calculate Q, we first note that for each source 
with a given luminosity, L, and colour, C, there is a redshift, 
^max, at which the observed flux is equal to the detection 
threshold 



S{L, C, Ao, 2max) — Sth — 



(l+2max)I/^(I/,C,Ao) 
^TvDliZn,,^) 



(14) 



where Dl is the luminosity distance and L^ is the rest- 



Now that we have all the tools ready, we can test the 
capability of our model in reproducing the observed prop- 
erties of 850 ^m sources and find which particular choices 
of model parameters are implied. In the following sections, 
we first discuss the observational constraints that we want 
to reproduce with our model and then we will proceed with 
finding the best-fit model which can reproduce those prop- 
erties. 



3 850 ^M OBSERVATIONAL CONSTRAINTS 

3.1 Observed 850 /im source count 

Among several existing extragalactic s ubmm surveys which 
provide source counts at 850 /im ( Coppin et al.l l2006l : 
IWeiss et all l2009l : lAustermann et all I2OIGI ). the SCUBA 
Half- Degree Extragalactic Survey (SHADES) (jCoppin et al.1 
12009) is the largest one which has the most complete and un- 
biased sample. However, this survey and other surveys which 
use JCMT and the same blank field method are restricted 
by the JCMT confusion limit of ~ 2mJy at 850 /xm and 
cannot probe the source counts of the fainter population. 
Using a complementary method, the lensing technique has 
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Figure 2. A compilation of some observed 850 /^m source counts. 
Blue circles, gree n triangles and red squares are respectively data 
points taken from Knudsen et al.l | |200S|) . ICoppin et al.l 1120061 ) and 
Ijohansson et al.l l|201ll 'l. 



been used to probe 850 /xm source co u nts to flux thresholds 
as low as O.lmJv (ISmail et al.1 12002| : iKnudsen et"ai] l2008l : 
[Johansson et al.ll201ll ). As we will show later in this section, 
the sensitivity of bright and faint submm source counts to 
the evolution of dusty galaxies is completely different and 
it is essential to incorporate a large dynamic range to con- 
strain possible evolutionary scenarios. Therefore, we use the 
best available observational information at both faint and 
bright tails of 850 /im source count, by combining al l of the 
SHADES data points with those of iKnudsen et al.l (J2008l l 
at flux thresholds < 2mJy to assemble our reference source 
count which is also in agreement with other observations 
(see Figure [2]). 



3.2 Redshift distribution of bright 850 /im sources 

Although the observed number counts at 850 fim, especially 
at faint fluxes, are sensitive to the redshift distribution of in- 
frared galaxies, they do not constrain it directly. As we will 
discuss later, different models with different redshift distri- 
butions can reproduce the 850 /im source counts with the 
same accuracy. Therefore, it is important to impose an ad- 
ditional constraint on redshift distribution of those objects. 
Unfortunately, studying the redshift distribution of 
submm galaxies is extremely difficult and there are only a 
few works on spectroscopically confirmed redshift distribu- 
tion of bright subm m galaxies with observed fluxes grater 
than ~ 4 — 5 mJy (jChapman et alJ l2005l : IWardlow et al.l 
I2OIII ). Those studies show a redshift distribution peaking 
around 2 ~ 2 (see the right panel in Figure [6|. We use this 
redshift distribution as one of our observed constraints and 
force our model to reproduce a redshift distribution which 
peaks at the same redshift. 



4 FINDING 850 mM BEST-FIT MODEL 

In order to be able to extract meaningful trends and differen- 
tiate between them, it is important to keep our phenomeno- 
logical model as simple as possible. It is therefore desirable 
to find the minimum number of free parameters without 
which the model cannot produce an acceptable result. Fur- 
thermore, the role of different parameters are not identical: 
while some parameters are not strongly constrained, small 
changes in others can change the result significantly. It is 
also important to look at degeneracies between various pa- 
rameters; some parameters are not independent and varying 
one may be compensated by varying the others. In this sec- 
tion we investigate our model to understand those issues 
before presenting our best-fit model. 



4.1 The source count curve: Amplitude vs. Shape 

To match the observed source counts, our model should be 
able to reproduce both the typical number of sources (i.e., 
the amplitude of the source count curve as a function of 
fiux threshold) and the ratio between the source counts at 
faint and bright fiux thresholds (i.e., the shape of the source 
counts curve). Following this line of argument, we can cate- 
gorize our model parameters into two different groups: those 
which play a stronger role in forming the amplitude of the 
curve and those which mainly affect its shape. As mentioned 
earlier, the luminosity evolution (see equation [T]) has the 
dominant role in determining the amplitude of the source 
count curve while its shape is controlled by other ingredients 
of our model which are colour evolution and the evolution 
of LF slopes. 

These trends are illustrated in Figure [S] In the left 
panel, the relative contribution of each parameter in chang- 
ing the amplitude of the source count curve is visualized; the 
length of the coloured segments along each axis (coloured 
segments along different axis are connected to each other) 
represents the sensitivity of source count amplitude to that 
parameter. The length of coloured segments is computed by 
changing all the model parameters (one at a time) by the 
same fraction, for instance 10%, and measure the change 
caused in the total source count. The segments which are 
connected by dot-dashed (red) and solid (blue) lines corre- 
spond respectively to a decreasing and an increasing param- 
eter. In the right panel, the relative contribution of differ- 
ent parameters in controlling the shape of the source counts 
curve is shown. In this figure, the length of coloured seg- 
ments represents the change in the ratio between faint and 
bright source counts. This time, we measured how much the 
ratio between a very faint source count, like 1 mJy, and a 
very bright one, like 1 Jy, is changing due to a fixed change 
in different parameters (e.g., 10%). The dot-dashed (red) 
and solid (blue) lines connect the coloured segments which 
correspond to an increase or a decrease in the values of the 
parameters. 

As is evident from these diagrams, the amplitude of the 
source count curve is strongly sensitive to the parameters 
of the luminosity evolution, particularly Za and n, while its 
shape is determined mainly by an evolution in the shape of 
LF and/or the strength of the colour evolution. 



© 0000 RAS, MNRAS 000, 000-000 



Genesis of the dusty Universe: modeling suhmillimetre source counts 7 







Amplitude 




z - 

b 


1.Z . z 
\ a / °' 


a 




\/\"x „ 


a.^ 




/\ 




W 


m 







Shape 








z — z 

6 \ a 


J 


f\ 




\,^ 


^^^^^* 


\/\\ 




n 


^^^^^ 


/\ /' 








/ 


rX 








/ 

It; 




m 





Figure 3. The sensitivity of the model to various parameters. While the left panel illustrates how much the variation of different 
parameters could change the amplitude of the 850 /im source count curve, the right panel shows their effect in changing the shape of the 
source count curve. To measure the amplitude changes, we varied each parameter by 10% around its best-fit value and measured the 
difference in the total source counts. The shape sensitivity is probed by measuring the relative change each varying parameter can cause 



in the ratio between two source counts at faint and bright flux thresholds (i.e. 



^lm.Iy 
JVlJy 



). The lengths of the coloured segments on each 



axis shows how large those changes are. The blue (red) segments which are connected by solid (dash-dotted) lines represent the result 
as we increased (decreased) the value of each parameter. We show parameters a^ and a^g on one axis since their effect in changing the 
source count properties is similar. 







Figure 4. The likelihood (i.e. exp{—^x^)) distribution of different parameters for the best-fit model. For each panel, we fixed all the 
parameters at their best-fit values and changed only one parameter at a time. By varying each parameter, the quality of the fit and 
hence the likelihood is changing. We normalized the amplitude of curves such that the area under each curve (i.e. the total probability) 
be unity. 



4.2 The luminosity evolution 

The luminosity evolution is the backbone of our model and 
has the most important role in producing the observed 
source counts and the redshift distribution of 850 ^m ob- 
jects. However, it is not surprising that the calculated source 



count is not strongly sensitive to the model properties at 
high redshifts. This is mainly because objects at very high 
redshifts have decreasing fluxes (despite the advantageous 
ii'-correction). Consequently, the chosen value of m has a 
negligible effect on the amplitude and shape of the source 
count curve (see Figure [Sjl . Nevertheless, one should note 
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be compensated by a very small change in n and/or Za- 




Figure 5. Map of likelihoods for models with different n and Za 
where [w,aa,ap] = [0,0,0]. The region which contains the best- 
fit models is in the middle of the coloured band and surrounded 
by the dashed blue contours which indicate the likelihood of 0.85. 



that this argument would not necessarily work if the popu- 
lation of bright IR galaxies continued to evolve with look- 
back time for all redshifts. In other words, the exact value of 
m is not strongly constrained by the observed source counts 
and all the negative values can produce similar results (see 
the middle panel in the bottom row of Figure |4]) . 

Contrary to m, the source count curve is very sensitive 
to the exact values of n and Za which respectively control 
the rate by which the characteristic luminosity of IR galaxies 
increases with redshift and up to which redshift this growth 
continues. However, there are two solutions for reproducing 
the same source counts: one is to increase the characteristic 
IR luminosity rapidly up to a relatively low redshift and the 
other one is to increase the characteristic IR luminosity by 
a moderate rate but for a longer period of time (i.e., up to 
higher redshifts). Therefore, as illustrated in Figure [S] there 
is a degeneracy between n and Za and one cannot constrain 
them individually by looking at the observed source counts. 
However, additional information about the redshift distri- 
bution of submm galaxies or the slope by which the char- 
acteristic IR luminosity is growing at low redshifts could 
resolve this degeneracy. As we will discuss later, the former 
constraint is used in finding our best-fit model. 

Unlike n and Za, the length of the redshift interval 
during which the luminosity evolution remains constant 
before starting to decline, Zb — Za, is not an essential 
part of our model. We only introduced this feature to 
have a smoother transition between growing and declining 
characteristic IR luminosity and also redshift distribution 
of submm galaxies. As Figures [3] and 0] show, variations 
in the adopted value for zi, — Za do not change the source 
count significantly. Moreover, any change in its value could 



Based on these considerations and as we already men- 
tioned in Section [2. 2 1 we reduce the number of free parame- 
ters we use in our model by choosing m = —1 and zt^Za ~ 1. 
Furthermore, the observed redshift distribution of submm 
galaxies which peak s around 2 ~ 2 IjChapman et al.ll2005l : 
IWardlow et al.ll201ll ). limits the acceptable values of Za to 
~ 1.6 (see the right panel in Figure [BJ. 



4.3 Other necessary model ingredients 

Although the luminosity evolution is necessary for produc- 
ing correct number of observable sources, it is not sufficient 
to provide a correct shape for the source count curve. Based 
on our experiments by varying different parameters of the lu- 
minosity evolution, one can either get a good fit at the faint 
number counts and under-produce the bright end or pro- 
duce correctly the bright source counts and over-estimate 
the faint objects. In other words, other model ingredients 
like the colour evolution and/or the evolution of LF slopes 
are required to adjust the shape of the source count curve ap- 
propriately in order to fit the faint and bright source counts 
at the same time. For instance, the colour evolution can 
be used to help the model with under-production of bright 
sources and the evolution of LF slopes can compensate for 
the over-production of the faint sources. 

As a starting point, we have first explored fits using ei- 
ther the colour evolution or the evolution of the LF slopes, 
in order to keep the model as simple as possible. We did this 
by trying to fit the 850 fim source counts once using a combi- 
nation of the luminosity evolution together with the colour 
evolution when the slopes of LF do not evolve (i.e., No- 
a-Evol) and once using the luminosity evolution combined 
with the evolving LF slopes, without evolving the colour dis- 
tribution (i.e. No-C-Evol). The results are shown in FigurelH] 
and are compared with the best-fit result when the colour 
evolution and the evolution of LF slopes are both active 
(i.e., 850-model). The quality of the fit for "850-model" is 
better than both of the other cases which may be attributed 
to the additional free parameter used in the "850-model". 
However, both "No-a-Evol" and "No-C-Evol" models have 
unfavourable implications. The best-fit for the "No-C-Evol" 
case (see Table [1} requires a large increasing slope in the 
luminosity evolution function which gives rise to a viola- 
tion of the counts of IR sources at short wavelengths (e.g. 
at 70 /im, see Section [5} . The "No-a-Evol" fit on the other 
hand, requires a steep colour evolution which is too extreme 
to be acceptable, since it would imply that at redshifts 2 > 1 
all infrared sources, independent of their luminosities, have 
the same dust temperatures as low as T ~ 20K. However, 
when both colour evolution and evolution of the bright and 
faint-end slopes of the LF are allowed simultaneously, these 
problems disappear which makes this solution preferable. 



4.4 The 850 nm best-fit model 

As we showed, it is necessary to have an evolutionary model 
which incorporates a luminosity evolution, colour evolution 
and also evolving LF slopes, in order to fit the observed 
850 ^m source counts. We also argued that we can fix some 



© 0000 RAS, MNRAS 000, 000-000 



Genesis of the dusty Universe: modeling suhmillimetre source counts 9 



10' 



10" 



10' 



10^ 



10' 





850-modeL 


~. 


..^ 


No-C-Evol ■ 






f]^>^^ No-a-Evol 






^"y^-i^ot- ^ 


No-z.-Const 




''■■■■T^^^ 








" 


"'^^^ 














_ 


^^. _ 




^^vy 




^V* 




"^ %i 




3^ 




■^T^ 


.. 


- 


% Knudsen et at. 2008 '^ji^ _| 




w Cop-pin et aL 2006 












-L 






. . .>->s.l 



0.4 



0.3 



0.2 



0.1 



0.0 



Chapman et at. 2005 



850— model 

l\U-L. —HjUOL 

No — a—Evol 
No-z^-Const 



1.0 



10.0 



100.0 



S (Jy) 




Figure 6. Left panel: The best-fit model which is constrained by 850 ^m source count and the redshift distribution of subnim galaxies 
(the blue solid line) is illustrated next to the reference observed data points (see Section [S}. For comparison, three other best-fit models 
are also plotted: the best-fit model without colour evolution shown by the dashed orange line (i.e., "No-C-Evol"), a model without 
evolving LF slopes shown by the dot-dashed green curve (i.e., "No-a-Evol") and finally the best-fit model which is only constrained by 
the source count is shown using the dott ed blue curve. Right p anel: the comparison between the observed redshift distribution of 850 fim 
galaxies which are brighter than 5mjy llChapman et al.ll2005l ) and what our best-fit model implies. The histogram shown in red with 
dashed line is the observed probability distribution and the histogram with solid blue line shows the probability distribution of similar 
objects in our model. For comparison, the best-fit model without constrained redshift distribution, is shown with dotted blue histogram 
("No-Za-Const"). The other two models, "No-C-Evol" and "No-a-Evol" are respectively shown by dashed orange and dot-dashed green 
lines. 



Model 


n 


Za 


w 


a.a 


a/3 


850-model 


2.0 


1.6 


2.0 


0.6 


0.4 


No-C-Evol 


2.8 


1.6 





0.6 


0.6 


No-a-Evol 


1.6 


1.6 


5.6 








No-^a-Const 


2.4 


3.6 


2.4 


1.0 


2.2 



Table 1. Parameters which define different best-fit models constrained to reproduce the observed source counts at 850 /^m. All the 
models are using Z;, — Za = 1 and m = — 1 and except the "No-Za-Const" model, all of them are constrained to reproduce the redshift 
distribution of submm galaxies and therefore use Za = 1.6. The predicted source count each of those models and their implied redshift 
distribution for submm sources is illustrated in Figure El 



of the initial parameters of the model since they have no 
significant effect on the results and chose proper values for 
them (i.e., m = — 1 and Zb — Za = 1) ■ Moreover, we showed 
that we need to set the redshift at which the luminosity 
evolution peaks to Za = 1.6 to reproduce the redshift distri- 
bution of submm galaxies correctly which also resolves the 
degeneracy between Za and n. After taking into account all 
of those considerations, we end up with 4 free parameters in 
our model which are needed to be adjusted properly to re- 
produce our observed sample of 850 /^m source counts; those 
parameters are the slope of the colour evolution, w, the rate 
by which the faint and bright end slope of LF is changing 
with redshift, a^ and a^, and finally the growth rate of the 
characteristic IR luminosity, n. 

Since our algorithm for calculating the source count 
is fast and accurate, contrary to much more cumbersome 
Monte-Carlo-based approaches, we can perform a compre- 



hensive search in the parameter space for the best-fit model 
instead of choosing it "by hand" . We split each dimension 
of relevant regions of parameter space into equally spaced 
grids and calculate the source count for modes associated 
with each node of our grid structure. Then we calculate the 
likelihood of each model for reproducing the observed source 
counts (i.e. exp{—^x^))- Finally, we choose the model with 
maximum likelihood as our best-fit model. The parameters 
which define the best-fit model for 850 /.im source counts, 
"850-model" , is shown in Table [l] together with those of 
"No-C-Evol" and "No-a-Evol" which we discussed in the 
previous section. Those models are also compared with ob- 
servational data sets used for constraining them, in Figure [6] 
where also the redshift distribution implied by each model 
is illustrated (in the right panel). 

It is also interesting to inspect the properties of a best- 
fit model in which the peak of the luminosity evolution, Za, 
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is not fixed. The parameters which define this model are 
presented in the last row of Table [T] and its source count 
and redshift distribution are shown by the blue dotted lines 
in Figure [6] The quality of the fit for this model is even 
better than "850-model" except at very low flux thresholds 
(unsurprisingly, since it has an additional free parameter, 
namely Za). The rising slope of luminosity evolution, n, and 
the colour evolution, w, in this model are not drastically 
different from the "850-model" but its peak of luminosity 
evolution happens at higher redshifts and the model requires 
much steeper slope evolutions to compensate for too many 
observable sources. Not only is this model unable to produce 
the observed redshift distribution ( see the right panel in 
Figure |6]), it predicts too few observable sources at shorter 
wavelengths which makes it unfavorable. 



5 OTHER WAVELENGTHS 

In the previous section, we discussed our model parame- 
ters and their role in producing the observed 850 /^m source 
counts and the redshift distribution of submm galaxies. We 
used those observed quantities as constraints to find our 
best-fit model, "850-model" . If we assume the evolution sce- 
nario that our best-fit model suggests is correct, and also 
the set of SEDs we used for reproducing the 850 /im source 
count are good representatives for real galaxies, then we ex- 
pect the same model, together with the same set of SEDs, 
to reproduce the observed number counts of IR sources at 
other wavelengths. Moreover, one can use the source count 
at other wavelengths to find a best-fit model for that specific 
wavelength and again expects to reproduce the source count 
at other wavelengths correctly. 

Before examining those expectaions, it is important to 
note that there are some fundamental differences between 
the source counts at shorter wavelengths and in the submm. 
Most importantly, as illustrated in Figures [7| and IIOI ob- 
jects observed at short wavelengths (e.g., 70 /im) are mainly 
at low redshifts while the observed submm galaxies are dis- 
tributed in a wider redshift interval and at higher typical 
redshifts. Also, the brightest objects at 70 /^m are the clos- 
est ones which is n ot the case for brightest submm galaxies 
(jBerta et al.ll2011f ). Moreover, there are additional physical 
processes, like AGNs, which can change the energy output 
of galaxies at shorter wavelengths. Since we do not take 
into account AGNs as a separate population, our model 
is not expected to reproduce necessarily good results for 
wavelengths shorter than 60 — 70 fim (see also the discussion 
in Section [2. 3p . Going from longer to shorter wavelengths, 
the observed sources are typically at lower and lower red- 
shifts. Consequently, the source counts at short wavelengths 
are only sensitive to the very low redshift properties of our 
model. In fact, the 70/im source count is mainly sensitive 
to the growth rate of the luminosity evolution, n, and the 
evolution in LF slopes; but by going to longer wavelengths, 
the colour evolution and the redshift at which the luminosity 
evolution is peaking, Za, become more important. In other 



^ In fact the wide distribution of observed submm galaxies in red- 
shift space, combined with the availability of a measured redshift 
distribution are the primary reasons we chose their observational 
properties to constrain our model. 



words, if a model reproduces the observed source counts at 
short wavelengths, this forms a confirmation that the evolu- 
tion at low redshifts is represented correctly but for a best-fit 
model constrained by observations at shorter wavelengths, 
the properties of model at intermediate and high redshifts 
are not strongly constrained. 

To investigate the performance of different best-fit mod- 
els constrained by source count observations of wavelengths 
other than 850 /im, we use the same procedure we incorpo- 
rated in finding the "850-moder' and only vary the effective 
parameters of the model, namely n, w, aa and a/3, to fit the 
observed data. The parameters which define best-fit models 
at different wavelengths are shown in Table [2] and the source 
counts they produce at different wavelengths are illustrated 
next to the observational data points in Figure [S] Different 
panels are for source count at different wavelengths rang- 
ing from 1100 /im on the top left to 70 /im on the bottom 
right. We do not show the results for 350 ^m to maintain 
the symmetry of figures since for this wavelength models 
and their comparison with observed data are in many re- 
spects identical to the case of 250 fim. In Figure[51 the source 
counts produced by different models are also shown using 
lines with different styles and colours: purple short-dashed 
for 70 /im, green long-dashed for 160 /im, orange dot-dot- 
dashed for 500 /im and finally solid blue lines for 850 /im (i.e 
the "850-moder'). In the following we discuss those results 
by categorizing them in different wavelength ranges, namely 
850 /im and 1100 /im as long submm wavelengths, 500 /im, 
350 /im and 250 /im as SPIRE or intermediate wavelengths 
and finally 70 /im and 160 /im as short wavelengths. 

5.1 Long submm wavelengths: 850 /im and 1100 /im 

As we mentioned earlier, the SED of star forming galax- 
ies at long submm ranges is essentially controlled by the 
Rayleigh- Jeans tail of the dust emission which simply falls 
off smoothly (see Figure [TJ. This means that if a model re- 
produces the observed counts at 850 /im, a good fit to the 
observed counts at similar and longer wavelengths is guar- 
anteed that given the distribution of sources which produce 
those counts have a similar redshift distributions, which in 
turn makes those counts equally sensitive to different pa- 
rameters in our model. As illustrated in the top panels of 
Figure [HI both models which are constrained by observed 
70 /iin and 850 /im counts and produce a good fit to 850 /im 
data, also produce a good fit to 1100 /im source counts. Our 
experiments with other models also confirm that the qual- 
ity of the fit they produce for observed counts at 850 /im 
and 1100 /im is highly correlated. Both of the models con- 
strained by 160 /im and SPIRE (e.g. 500 /im) source counts 
over produce the long submm wavelength source counts: the 
" 160- model" over-produces submm source counts mainly at 
bright fiux thresholds due to its extreme colour evolution 
(see Table[21) but fits the observations at faint fiuxes, thanks 
to its luminosity evolution which is similar to those of " 850- 
model" and "70-model". However, SPIRE constrained mod- 
els (e.g. "500- model") have stronger evolutions both in terms 
of luminosity and colour and over-produce the data in all 
observed fluxes. 

One should note that t he 1100 /im data points are highly 
incomplete below 3mJy (JHatsukade et al.l 120111 ) and are 
prone to large field to field variations and relatively large 
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Figure 7. The modeled redshift distribution of objects which are observed with different flux thresholds at 70 /jm (in the left panel) 
and 850 fim (in the right panel). Three different flux thresholds which are 0.1, 1 and lOmJy are shown respectively using solid{blue), 
dot-dashed (green) and dashed (red) lines. 



A 


n 


w 


aa 


a/3 


850 Aim 


2.0 


2.0 


0.6 


0.4 


500 Atm 


2.6 


3.8 


0.0 


0.6 


350 Atm 


3.0 


4.4 


-0.2 


1.4 


250 Aim 


3.0 


4.4 


0.0 


1.4 


160 Aim 


2.2 


4.2 


0.8 


0.2 


70 Aim 


2.2 


2.0 


0.6 


0.6 



Table 2. Parameters which define different best-fit models constrained to reproduce the observed source counts at different wavelengths. 
The first column, A, indicates the wavelength for which the model is constrained to produce the best fit to the observed source counts. 
All the models are forced to reproduce the redshift distribution of submm galaxies and therefore use z^ = 1-6, Zj — z^ = 1 and m = — 1. 
The predicted source count each of those models is implying for different wavelengths is illustrated in Figure [8] 



errors in the whole range of observed fluxes. Moreover, the 
observed counts are at bright flux thresholds which makes 
them sensitive to a narrower redshift interval in comparison 
to the fainter flux thresholds (see the right panel of Fig- 
ure [7|) and because of those reasons they cannot constrain 
models better than what 850Atni source counts are capable 
of. Therefore, we only consider the "850- model" as a model 
constrained by long submm source counts. 



5.2 SPIRE intermediate wavelengths: SOOAini, 
350 Aini and 250 A^ni 

The three 500 Ainr, 350 Aim and 250 Ainr wavelengths are close 
together and besides being in the middle of wavelength 
ranges we study, they have intermediate properties with re- 
spect to the redshift range each wavelength is mostly sensi- 
tive to: as it is shown in Figure [TOl while at the bright flux 
thresholds the source counts mainly consist of low redshift 
sources, at fainter fluxes they are sensitive to the interme- 
diate (i.e. 1 < 2 < 2) and high redshifts (i.e. 2 < z). Natu- 
rally the general behavior of models for 500 nm is closer to 
850 Aim while 250 fim. is close to shorter wavelengths. How- 
ever, those intermediate wavelengths are closely similar to 



each other more than being similar to other wavelengths. 
Consequently, the best-flt models constrained by SPIRE 
wavelengths have similar parameters ( 250-model and 350- 
model have almost identical parameters) and any of those 
models agrees with the observed source counts of the other 
two. However, SPIRE-constrained models all require steep 
luminosity evolution and strong colour evolution in addi- 
tion to a steepening bright end slope of the LF with red- 
shift (i.e. positive Qq together with almost zero ap, see Ta- 
ble [2]). The first two mechanisms produce too many objects 
in comparison to what is needed for the source counts at 
other wavelengths while the steeper bright end of the LF 
partially compensate the over-production of bright objects: 
as is shown in Figure [S] SPIRE-constrained models over- 
produce the faint source counts both at longer and shorter 
wavelengths but the steep bright end of LF produces results 
which are close to the observed faint source counts at very 
long wavelengths (i.e. 850Aini and 1100 Aim) which in turn 
is too strong to leave enough sources required for the bright 
160 nm sources. As we will discuss later, it is not surprising 
that SPIRE- models agree with bright 70 fj.ra counts since 
they are essentially z ~ objects for which evolutionary 
mechanisms in the model barely have any effect. 
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Figure 8. In different panels, the best-fit model predictions for cumulative source counts at 850 ^m and differenti al source counts at 1100, 
500, 250, 1 60 and 70 /jm are plotted next to the observed data. The 1100 ^m observational d ata sets are fromlScott et al.l 1120101 ) (red 
diamonds) lAustermann et al.l 1120 091. l2010h ( orange circles and green squares respeci tvely) and iHatsuka de et a l.l 11201 If) (blue t riangles) ; 
data at 850 /xm is from I Coppin et al.l 1 120061) (green triangles). iKnudsen et al.l 1120081 ) (red circl es) and [Johansson et al.l 1 120111) (orange 
squares) . The data for 500 and 250 ^tm source counts are from the B LAST experiment take n from lPatanchon et al.l l(2009|) (red diamonds) , 
Herschel data taken from lOIiver et al.l 1120101 ) (green squares) and [ Clements et al.l 1120101 ) (orange circles). We also used one data point 
at 350 /im from the SHARC 2 survay (green triangle) JKhan et al.ll2007l ). The data at 160 and 70 /xm are based on Spitzer observations 
taken from lBethermin et al.l 1(2010) where the data points shown by filled squares represent stacking results. Different lines correspond 
to the source counts produced by various models which are constrained to fit the observed source counts at different wavelengths: blue 
solid line shows the "850- mode", purple dashed line shows the "70- model" and green long-dashed and orange dot-dot-dashed lines are 
for "160- model" and "500- model" respectively (see Table [2]l 
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While the SPIRE models fail to agree with observations 
at shorter and longer wavelengths, the two models which 
are successful at long submm ranges (i.e. the "850- model" 
and " 70- model" ) also have a reasonable agreement with the 
SPIRE data. However, they under-produce the counts at 
Sth < lOOmJy by a factor of ~ 2. As we will show later, as- 
suming the "SSO-model" to be correct, this under-production 
is a hint for a population of cold luminous IR galaxies resid- 
ing only in low and intermediate redshifts (i.e. z ~ 1). This 
is also consistent with the steep coluor evolution which is 
implied by best-fit models at those wavelengths which pro- 
duces enough cold galaxies at lower redshifts. 



5.3 Short wavelengths: 160 /im and 70 fim 

At short wavelengths like 70 /im, the K-correction is not 
strong enough to counteract the effect of cosmological dim- 
ming which makes the source counts at these wavelengths 
almost insensitive to the properties of IR galaxies at high 
redshift. Since the starting point of our model is the ob- 
served distribution of IRAS galaxies which are selected to 
be local galaxies with Seofim > Uy, and the SED templates 
we use are extensively tested to match these galaxies, any 
variation of our model by construction cannot disagree with 
bright 70 /im counts q However, at fainter flux thresholds 
(e.g. Sth < lOOmJy) the sensitivity of 70 /im counts to low 
and intermediate redshifts, 0.5 < z < 2, increases (see Fig- 
ure nop which makes them more sensitive to the slope of 
the luminosity evolution. As a result, all the models agree 
with 70 /im observations at bright flux threshold while the 
SPIRE-model which has steeper luminosity evolution (i.e. 
greater n, see Table ^ diverges from observed counts and 
over-produces them. It should be noted that the " 70-moder' 
which is tuned to agree with mainly low-z IR sources per- 
forms as good as "850- model" in longer wavelengths. 

At 160 /im, the redshift distribution of sources is almost 
identical to the distribution of galaxies which are responsi- 
ble for 70 /im: the bright end of the counts is governed by 
local IR galaxies while for flux thresholds Sth < lOOmJy the 
intermediate redshifts (i.e. 0.5 < z < 2) become important. 
However, the flux threshold at which the turn over between 
domination of local and farther sources happens is one order 
of magnitude higher for 160 /im sources (see bottom row in 
Figure fTH)) . Moreover, at 160 /im the domination of interme- 
diate redshift sources in shaping the faint source counts is 
slightly stronger than at 70 /im. In other words, the 160 /im 
source counts are slightly more sensitive to the distribution 
and evolution of distant IR galaxies in comparison to 70 /im 
counts. 

The best-fit model which is constrained by 160 /im 
source counts has similar slope of luminosity evolution, n, 
compared with the "850-moder' and "70-moder' but re- 
quires much stronger colour evolution and larger fraction of 
bright objects (i.e. steeper faint end slope and flatter bright 
end slope of LF at higher redshifts). 



The SPIRE-models on the other hand, violate the ob- 
served 160 /im counts by under-producing them at bright 
end and over-producing them at faint fluxes. The "850- 
moder' and "70- model" on the other hand, match the 
faintest 160 /im count but under produce the brighter ob- 
jects by a factor of ~ 2 which is the same factor show- 
ing up in the difference between what those two mod- 
els produce and observed faint SPIRE counts. A similar 
discrepancy between models and 1 60 /im counts has been 
pointed out in som e recent works (|Le Borgne et al, l2009l : 
IValiante et"aLll2009l '): while the lLe Borgne et al.1 (120091 ') best- 
fit model which produces correct source counts at 70 /im and 
850 /im deviates from observati ons between lOmJy < Sieo < 
lOOmJy by a factor of ~ 2, the lVaUante et al.l (|2009l ') model 
underestimates number counts at 160 /im by a factor of ~ 5. 
We discuss this issue further inl5.4l 



5.4 A best-fit model for all wavelengths 

As we showed in previous subsections, the required evolution 
scenarios for different models which are constrained by var- 
ious wavelengths is too diverse to be reconciled in a single 
model; however, the 70 /im and 850 /im source counts can 
be explained by the same model, either the "850- model" 
or "70-model" which are almost identical in terms of im- 
plied evolution and distribution of IR galaxies and predicted 
source counts at different wavelengths. While 70 /im source 
counts are produced mainly by local and low redshift ob- 
jects, 850 /im sources are distributed in higher redshifts and 
in a wide redshift interval. This means the model which can 
explain the source count at both 70 /im and 850 /im gives a 
consistent distribution of IR galaxies from local Universe to 
very high redshifts. In other words, the same evolution sce- 
nario for star-forming galaxies which we observe locally, can 
also produce the observed 850 /im counts. This scenario, con- 
trary to what is implied by the best-fit models constrained 
at other intermediate wavelengths, does not strongly vio- 
late observed source counts at other bands, though it under- 
produces them at some flux thresholds. On the other hand, 
models which agree with observations at intermediate wave- 
lengths require extreme color evolutions to produce larger 
number of cold objects in intermediate redshift ranges which 
is relevant for the source count at those wavelengths (i.e. 
1 < z < 2, see the left panel in the middle of Figure fUJ)) . 
However, this steep colour evolution at higher redshifts (i.e. 
z > 2) produces too many observable sources at 850 /im. 

Moreover, as we mentioned earlier different authors 
with fundamentally different models have reported in- 
consistencies in their best-fit models (which fit the long 
and short wav elengths) and interme diate source counts 
(|Le Borgne et a l. 2009; Valiante et all 2009 ) Fl and the un- 
derestimation of 160 /im (and SPIRE wavelengths) seems to 
be model independent. 

There are two main possibilities which can explain the 
origin of discrepancy between the "850-model" and source 
counts at SPIRE band and 160 /im, besides doubting the va- 
lidity of observed counts: (i) the IR SEDs in our model are 



^ Unless the luminosity of IR galaxies increase with redshift ex- 
tremely steep to produce high-z objects which are observed in 
70 /im band as bright as the brightest local IR galaxies despite 
the cosmological dimming and negative K-correction. 



•^ Those works point out this issue only for the source counts at 
160 /tm since the observations at 250 /im, 350 /tm and 500 /tm have 
been available only recently. 
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Figure 9. The same set of observed source counts which is shown in Figure|8]is also illustrated here. The blue solid line shows the result 
of "850- model" using a modified set of SED (see Secti on 15.411 to correct un der-p roduction of bright 1 60 /.tm sources and faint SPIRE 
source. Source count produced by two different mode ls. iLagache et al.l ll2004h and IValiante et al.l l|2009f l , are also respectively shown by 
orange dash ed line and gray dot- dashed lines. While I Valiante et al.l 1120091 ) and "850-moder' counts at 160 /^m are almost identical, we 
did not find lLagache et al.l l|2004l 'l model predictions at 1100 /.tm to include them in the top-left panel. 
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Figure 10. The fractional contribution of sources in different redsfiift bins in producing the observed source counts at different wave- 
lengths are illustrated. The blue Dotted, cyan dot-dashed, orange dot-dot-dashed and red dashed lines respectively show the source count 
produced by objects in < 2: < 0.5, 0.5 <z<l, 1<2<2 and 2 < z < 5 redshift intervals while the total source count is shown using 
the brown solid line. In the top section of each panel, the fraction of sources for each redshift bin in producing source counts at different 
flux thresholds is plotted. 
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not representative and should be modified in a certain way 
to produce more flux at observed wavelengths Xots ~ 160 
- 500 fim or (ii) existence of a population of cold galax- 
ies at relatively low redshifts. Indeed, there is a body of 
evidences supporting a population of cold galaxies which 
are underrepres ented in IRAS gal a xies a r id hence 70 ^tm 
source counts llStickel et al. 1998I. 2000l: Chapman et al 



_2002; P atris et ai.l '2003':'Dcnnefe ld et al.ll2005l : ISaiina et al 
2006 : Amblard et al. 2010 ). However, based on available 



data there is no possible way to disentang l e betw een the two 
mentioned possibilities (ILe Borgne et al.l 120091) . Therefore, 
due to its simplicity, we try to find a modification in SED 
amplitudes which could improve the agreement between the 
"850-model" and observed source counts at SPIRE wave- 
lengths and 160 ^m, without aflecting other wavelengths. In 
the following we first introduce the desired SED modifica- 
tion and then based on the redshift distribution of sources 
responsible for different source counts in our modified model, 
we try to constrain different properties a population of cold 
IR galaxies should possess to be equivalent to that SED 
modification. 

Since the "850-model" is capable of fitting the observed 
70 /xm source counts, any modification of SED templates 
should be at wavelengths longer than 70 /im to leave this 
agreement intact. The 850 fira source counts that model pro- 
duces should also remains the same which put an upper 
limit for the allowed wavelength range of any modification: 
since a significant fraction of observed counts at 850 ^im 
is produced by sources at high redshifts, typically z > 2 
(see Figures [7] and [5| , SED templates should not change at 
rest-frame wavelengths longer than ~ 200 /im. Based on this 
ansatz, we searched for the amplitude and the wavelength 
range of SED changes which can bring the "850- model" in- 
termediate wavelength counts close to the observed values. 
We found that if the SED templates we use in our model 
be amplified by a factor of 1.6 in the rest-frame range of 
70 < A < 150 /im, the "850-moder' can fit the SPIRE and 
160 /im source counts, without any change in already de- 
cent results at 70pim, 850 /xm and 1100 /xm. Figure O illus- 
trates the results produced by this modified "850-moder' 
at different wavelengths (b l ue so lid line) together with the 
results from J Lagache et al.l (|200J) (orange dashed line) and 
IValiante et al.l 1 20091 ) (gray dash-doted line) models. The 
redshift distribution of sources responsible for those results 
are also shown in Figure [TIJl 

The SED boost we implemented, recovers the observed 
160 nm counts by doubling the number of observable sources 
with 5160 > 10m Jy; as illustrated in Figure [^ a signif- 
icant fraction of those sources is distributed at redshifts 
z < 1 which is a l so in agreement with recent observa- 
tions ({Jacobs et al.ll201ll : iBerta et al.l 120111 ). However, for 
the SPIRE source counts, the SED modification increases 
the number of observed sources with Sspire < lOOmJy 
which are mainly at redshifts 0.5 < z < 2. Those redshift 
distributions together with the fact that the emission from 
cold dust with temperatures T < 30K peaks at ~ 100 /im, 
make our SED modification equivalent to adding a popula- 
tion of preferentially cold dusty galaxies which do not ex- 
ist at high redshifts. The width of redshift range in which 
those objects can exist depends on their temperature range: 
a warmer population can have a wider redshift distribution 
but should have a typically higher redshifts to be invisible 



at 70 /im while a colder population can only exist in low 
redshift ranges in order to not interfere with 850 /im counts. 
Interestingly, this is in agreement with the required steep 
colour evolution implied by models constrained by 160 /tm 
and SPIRE band counts (see Table [2]). 

As one c an se e in Figure [^ some m odels like 
iLagache et"aLl (|2004l ) and IVahante et~aLl l|2009l ) not only 
produce enough observable sources in SPIRE and 160 /im 
bands, but at some flux thresholds produce too many 
sources. However, we notice that both of those models 
need an additional mechanism to compensate for under- 
production of visible sources at intermediate wavelengths 
which mimic a po pulation of cold sources only in low red- 
shifts. For instance. iLagache et al.l (|2004l ) use a class of "cold 
galaxies" in their evolut ionary model which a re present only 
at low redshift, z < 0.5: iLag ache et al.l l|2003h show that this 
cold population is producing up to ~ 50% of the observed 
170 /iin flux which means without them the source count 
is reduced by the sa me factor, consistent with our flnding 
and a lso other models |Le Borgne et al.ll2009l : IValiante et al.l 
I2OO9I ). Similarly. I Valiante et al.l 1 20091) need to strongly mod- 
ify the colour distribution of low redshift galaxies to correct 
for under-producing 160 /im sources by a factor of ~ 5; they 
modifled the observed colour distribution of IRAS galaxies 
which they use as starting point (similar to our model, see 
Section I2.1|l only for low redshift (i.e. z < 1) objects, by 
assuming an asymmetric Gaussian distribution which is 7 
times broader on the "cold" side of distribut ion in compar- 
ison t o the "warm" side (see Section 4.4 in IValiante et al.l 
(120091 )): even with this extreme modification their model 
does not completely match the observed 160 /im counts in 
addition to under-producing luminous 70 /im sources. 

Finally, it is worth mentioning that our SED modifi- 
cation is not expected to change the best-fit models which 
are based on 70 /tm and 850 /im since we used those source 
counts as constraints for our SED change search. However we 
double checked this issue by using the modified SED set and 
repeating the search in parameter space for a best-fit model 
which is constrained by source counts at 70 /im, 160 /im, 
SPIRE bands and 850 /im and recovering the parameters 
which define the "850-moder', for the best-fit model. 



6 DISCUSSION 

In this section, we discuss different properties of our best-flt 
model. First we discuss different implications of our model 
for the evolving properties of the IR galaxies like their colour 
and luminosity distributions. Then, we compare our model 
with other existing models for IR and submm source counts. 



6.1 The implied evolution scenario for dusty 
galaxies 

Our best-fit model mimics the number density evolution of 
IR galaxies by employing a luminosity evolution together 
with changing slopes of LF with redshift. While the for- 
mer changes the amplitude of LF, the latter acts to change 
the shape of LF properly to reproduce a number density 
history which matches the observed FIR and submm source 
counts. Although the good agreement between source counts 
which our best-fit model produces at different wavelengths 
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Figure 11. Top: The evolution of CF (left) and LF (right) with rcdshift. The lines from thick to thin are respectively for redshifts 2 = 
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and the observed numbers fLrmly supports our model, the 
implied results should be treated carefully, mainly due to 
the simplicity of the model: for instance, in our model we 
adopted the local LF of IR galaxies which is a simple dual 
power law function, and evolved its characteristic luminosity 
and slopes with redshift to distribute objects with different 
luminosities in redshift space correctly. 

Although the functional forms we chose for evolution 
of those parameters are rather arbitrary, the actual dis- 
tribution they produce at different redshifts is necessary 
to explain the properties of the observed sources. To 
emphasize this point, in Figure [11] we illustrate how the 
fractions of galaxies in different luminosity bins are evolving 
with redshift, together with the exact shape of the LF our 
model requires at different redshifts. As can be seen in the 
bottom right panel of Figure llll at low redshifts fainter 
objects dominate the population of IR galaxies but at 
higher redshifts (i.e., z > 2), objects which are in brighter 



luminosity ranges take over and become more dominant. 
Specifically, this diagram shows that beyond 2 ~ 1, objects 
in the Luminous Infrared Galaxies (LIRGs) class (with 
IR luminosities between 10^^ and 10^^ Lq) dominate the 
obscured cosmic energy production (although in terms of 
numbers, fainter galaxies still dominate). The diagram also 
shows that UltraLuminous Infrared Galaxies (ULIRGs, with 
IR luminosities exceeding 10^^ Lq), although increasing 
in importance towards high z, never dominate the cosmic 
energy production. These c onclusions are in agreement with 
analyse s of Spitzcr data by Le Floc'h et al.l (,2005i ) for z < 1 
and bv lMagneUi et all (120111 ) for z < 2.3. 

Our best-fit model implies a change in the shape of the 
LF with redshift and we are the first to consider this possibil- 
ity in a model of this type. This slope evolution implies that 
the faint end slope of the LF is flattening with increasing 
redshift (see the top right panel of Figure [TTT) . At first sight 
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it may seems surprising that our results are sensitive to the 
faint end slope at high redshifts, but the strong luminosity 
evolution combined with the strong negative ii'-correction 
brings sub-L, galaxies well within the region of detectabil- 
ity at 850 nm. Nevertheless, this result depends strongly on 
the sub-mjy counts at 850 fim which currently are derived 
from studies of gravitationally lensing clusters, and therefore 
sample limited cosmic volume. As such, these counts may 
be affected by cosmic variance and establishing their levels 
more firmly (e.g., with ALMA), is necessary for confirming 
our conclusion. However, we note the non-parametric esti- 
mates of the IR luminosity function at hig h redshifts (e.g., 
IChapman et al.ll2005l : IWardlow et al.ll201ll i do indicate flat- 
ter faint-end slopes than the local LF. These samples were 
selected at 850 ^m and may not be complete in IR luminos- 
ity, in particular they may be deficient in objects with high 
dust temperatures (e.g., Magdis ct al. 2010). In order to set- 
tle this point, high-z luminosity functions of IR luminosity- 
limited samples are required. Redshift surveys of Herschel- 
selected samples will be needed to construct such luminosity 
functions, and this will become feasible with facilities such 
as ALMA and CCAT. If confirmed, the flattening of the 
faint-end slope of the LF towards higher redshifts requires 
a physical explanation, which perhaps could be found using 
physically motivated models and simulations. One possibil- 
ity is that a higher metagalactic ultraviolet flux at higher 
redshifts would suppress the development of a star-forming 
interstellar medium in low-mass galaxies. Another possibil- 
ity is a stellar-mass-dependent evolution in Mdust / Msta.Ts to- 
wards higher redshifts. Such an evolution could result from 
the decreasing overall metallicity towards high z combined 
with the net effect over time of the buildup of dust as a result 
of stellar evolution and its consumption by star formation. 
The latter model can be tested observationally using a com- 
bination of Herschel data and multi-band optical imaging 
(Bourne et al., in prep.). 

We also showed that our model requires a colour evo- 
lution to reproduce the observed source counts. This colour 
evolution implies that objects with the same luminosities 
have lower typical dust temperatures at higher redshifts. In 
the bottom left panel of Figure [11] we showed how the frac- 
tion of objects with different temperatures is evolving with 
redshift for all the objects which have intrinsic luminosities 
between lO" < 7^ < lO". While at low redshifts a popu- 

lation of objects which have typical temperatureQ of 40K, 
dominates the colour distribution, at higher redshifts the 
cooler typical temperature of 30K is dominating. Moreover, 
while at low redshifts objects with warm dust temperatures 
are more common than the very cold objects, the situation 
is the opposite at earlier times. Although it is difficult to 
observe the evolution of dust towards cooler temperatures 
because of different selection effects, there is observational 
evidence for the implied col our (i.e., temperature) evolution 
llChapman et al.|[2 005: Pop e et al.ll2006 " 



1( 

lQ; 



Svmeonidis et al.l l2009i: iSevmour et al.l 
201(]| : lAmblard et al.ll201ol v 



Chapin et allbOOgl: 



2OI0I : iHwang et all 



The luminosity and colour evolutions discussed above, 

* We assign a single temperature to each SED based on its colour 
and finding a modified black body radiation with a fixed emissiv- 
ity, /3 = 1.5, which can produce that colour. 



show up in redshift distributions of source counts at differ- 
ent wavelengths. Figure [10] illustrates the buildup of source 
counts at different wavelengths by contribution from differ- 
ent redshift bins. It is evident that the longer wavelengths 
are reffecting the distribution of higher redshift dusty galax- 
ies while shorter wavelengths depend more on IR galaxies at 
lower redshifts. For instance 850 fj,Tn number counts mainly 
consist of galaxies with redshifts higher than 2: ~ 2, espe- 
cially at ffuxes ~ 10 — 20mJy; in addition, as shown in the 
bottom right panel of Figure 1111 at redshifts z ~ 2 the 
luminosity distribution of IR galaxies is not evolving con- 
siderably and even at higher redshifts the fraction of fainter 
sources, which have lower dust temperatures, goes up again 
in comparison to the brightest objects. Moreover, for two ob- 
jects with the same redshifts and intrinsic IR luminosities, 
the colder one will be brighter at observed 850 ^im. This 
means the visible 850 /xm sources at higher flux thresholds 
are preferentially colder than sources with fainter observed 
fluxes which are distributed typically at lower redshifts (see 
the increasing fraction of objects in 1 < z < 2 redshift bin 
going towards lower fluxes in the top right panel of Fig- 
ure llOf) . This is also in agreement with our experiments 
which show that at the bright end, 850 /im number counts 
are very sensitive to the slope of the colour evolution in our 
model, which is in fact the dominant mechanism to pro- 
duce those counts. Moreover, this implied broader tempera- 
ture distribution of observed sources at fainter flux thresh- 
olds and shorter wavelengths is in agreement with observa- 
tions which show bright 850 /xm population (i.e. > 4mJy) is 
highly biased towards cold dust temperatures while fainter 
sources (i.e. 1 — 4mJy) also contain an increasing fraction 
of more luminous objects with lower redshifts but warmer 
dust temperatures which make s them the main contributors 
to the 250 Atm source counts ijC hapman e t al.l |2004| . |201G| : 
iMagnelli et aLll2010l : [Casev et al.ll2 009, 2011). Consequently, 
at 2; ~ 2 the density of farlR selected ULIRGs is approxi- 
mately 2 times higher than that of 850 /^m-selected ULIRGs. 
This is consistent with our model: a typical ULIRG at red- 
shift 2 ~ 2 will be visible at 850 /im with a flux brighter 
than a few mjys and at 250 pim brighter than ~ lOOmJy; 
on the other hand, roughly ~ 20 — 30% of objects with few 
mjy fluxes at 850 /im are within redshifts 2 ~ 1 — 2 while 
around ^ 50% of objects with ~ lOOmJy fluxes at 250 /im 
are in the same redshift bin (see Figure [TD}. 

While the luminosity and colour distribution our best- 
fit model requires is consistent with observed 70 fim, 850 ^im 
and 1100 /im source counts, they are not sufficient to pro- 
duce enough sources at observed wavelengths in between. 
The inconsistency between models which successfully pro- 
duce 70 /im and 850 fim counts and their results at 160 /im, 
is also reported in other works and has been corrected 
by inc luding a popu lation of co l d galaxies at low red- 
shifts IJLagache et all 12001 |2004| : iLe Borgne et al.l l2009l : 
IValiante et al.l 120091 ). While it is important to note the 
under-production of 160 ^m counts in models, can be cor- 
rected equally by a modifying SED tem plates instead of 
introducing a new population (see also iLe Borgne et al.l 
l|2009r) ). there is some observational evidence for the ex- 
istence of a cold population at low redshifts which is 
under-represented in IRAS sa mple and is often associ- 
ated with bright spiral galaxies (|Stickel et al.l 119981 . I2OO0I : 
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Chapman et all I2OO2I: [Patris et al.1 l2003l: [Pennefeld et all 



2005l : ISaiina et al.ll20od : lAmblard et al.ll2010l ). 

Recently, flux density measurements at 250, 350 and 
500 /^m have become available for large samples of lo- 
cal galaxies from s urveys with SPIRE on the Herschel 
Space Observatory (|Eales et al.1 l2010l : lOliver et all l2010l : 
IClements et al.ll2010l V In addition to problems at 160 /xm, 
we also noticed the inconsistency between our best-fit model 
and the source counts provided by SPIRE observations. 
However, we resolved this issue by modifying our SED tem- 
plates to be able to reproduce the source count simultane- 
ously at 70 /im, 160 ^m, SPIRE band, 850 /im and 1100 /im. 
There is also observational evidence for a population of cold 
galaxies residing at low redshifts (equivalent to modified 
SEDs) in Herschel-selected samples: using a subsample with 
spectroscopic or reliable photometric redshifts from Herschel 
ATLAS survey, Amblard et al. (2010) performed isothermal 
graybody fits to low-redshift galaxies detected from 70 to 
500 fira, resulting in an IR luminosity-temperature relation 
offset to significantly lower tempera tures when compared t o 
the IRAS-based relation derived by Chapman et al.l ([2003). 
The new relation found by Amblard ct a l. (2010l) is consis- 
tent with earlier work bv lDve et al] (|2009l ') based on BLAST 
data. It is also important to realize that these results do 
not imply that the IRAS-based dust temperature fits are 
incorrect (in fact, they are often supplemented with mea- 
surements at longer wavelengths) but they imply that an 
IRAS-based selection is biased towards warmer objects. 

6.2 Our best-fit model and previous models 

There are several phenomenological models in the litera- 
ture which try to reproduce the properties of IR galaxies 
at differ ent wavelengths, with differ ent levels of complex - 



itv (e.g.. iBlain fc LongaiJ [19931: G uiderdoni et al.1 (^1997^: 



iBlain et a l. (1999): Ch arv fc Elbaj 120 01): Ro wan-Robinsoi] 



; [pole et al.,(,200l) : Eagache et al., (,2004'):lLewis et al.1 
LeBorgneerai] (|2009l ): IValiante et al.1 (|2009l )). In 



general, those models use an assumed form of luminosity 
evolution together with a density evolution to mimic the 
evolution of IR galaxy distributions. However, we have lim- 
ited ourselves to a pure luminosity evolution with no density 
evolution, and as pointed out in Section [2.21 the amount of 
density evolution allowed by the integrated CIB is small but 
does not have to be zero. On the other hand, we use evolving 
faint and bright ends slopes in our luminosity function. At 
first sight this may looks a simple substitution of free param- 
eters; we tried to investigate this by trying to substitute the 
slope evolution of our model with a density evolution. How- 
ever, the resulting fit with density evolution and in absence 
of slope evolution is significantly inferior to the fit obtained 
with pure luminosity evolution. 

Another usual practice in modeling the infrared and 
submm source counts is to use only one or a few SED tem- 
plates to represent the whole galajcy population. This ap- 
proach neglects the observed colour distribution of local IR 
galaxies and do not leave a ny possib i lity fo r colour evolution. 
However, similar to Valian te et al.l ([20091 ). in our model we 
use a complete set of SED templates which are shown to be 
representative of local IR galaxies. This choice enabled our 
model to explore the evolution in colours of IR galaxies in 
addition to their luminosities. 



Our model also differs from other works in the litera- 
ture in calculating the source counts for a given evolutionary 
scenario: we calculate the source count for a given model by 
computing the probability of observing different sources for 
different flux thresholds. The direct consequence of this new 
approach is a fast calculation routine which enables us to 
calculate the source counts at very bright observed fluxes, 
where Monte-Carlo based methods are very inefficient due 
to the rarity of such objects, and therefore sometimes end 
up with v ery noisy results (see the bright source counts pro- 
duced bv lValiante et al.l ([20091 ) in Figure (9]). Moreover, our 
fast algorithm is an important advantage when it comes to 
searching the parameter space for the best-fit model. 

A comparison between our best-fit model and 
iLagache et al.l ([20041 ). a s a model with o ut co lour distribu- 
tion and evolution and IValiante et al.l ([20091 ) as a model 
with colo ur distribution s is sh own in Figure (8] While at 
1100 /im, IValiante et al.l ([20091 ) produces ~ 2 times more 
visible objects than what our model produces, all models 
do reasonably well in accounting for the comula t ive 85 /im 
number counts. The results from IValiante et al.l ([20091 ) and 
our model are very close at SPIRE wavelengths but the 
iLagache et al.l (|200J) mo del overpredicts the b right counts 
at 250 /im. At 160 /xm, the lValiante et al.l ([20091 ) model over- 
produces the fain t counts and u nder- produces the bright 
objects. However, Lagache et al. (2004) fit the data better, 
while slightly over-produces the counts for intermediate to 
bright flux thresholds. Finally, at 70 /im, where all the mod- 
els are ex pected to fit the data s ince they use it as a starting 
point, the lValiante et al.l ([20091 ) source counts deviate from 
observations by over-producing the faint counts in expense 
of producing too few bright objects (probably because of 
a too extreme modification in colour distributions which is 
required in their model to compensate for a factor of ~ 5 
under-production of 160 /im sources). 



7 CONCLUSIONS 

We have described a backward evolution model for the IR 
galaxy population, with a small number of free parameters, 
emphasizing which parameters are constrained by which ob- 
servations. We also introduced a new algorithm for calculat- 
ing source counts for a given evolutionary model by direct 
integration of probability distributions which is faster than 
using Monte-Carlo sampling. This is an important advan- 
tage for searching large volumes of parameter space for the 
best-fit model. 

While most of the earlier works used only one or a hand- 
ful of SED templates to represent the whole population of 
IR objects, we used a library of IR SEDs which are able 
to match the IR properties of the large variety of observed 
star-forming objects. This approach is necessary in order 
to model the colour evolution of IR galaxies in addition to 
produce simultaneously the counts and the redshift distri- 
butions at wavelengths shorter than 850 /im. 

Contrary to some other models, we assumed a negli- 
gible contribution from AGN in our SED templates, not- 
ing the inclusion of AGN is only necessary for reproducing 
the properties of IR galaxies at very short IR wavelengthqj 

^ For instance at 24 /jm where our model under-produces the 
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which could also be sensitive to other modeling difficulties 
such as the PAH contribution to the SEDs. 

We used available 850 /^m source counts together with 
the redshift distribution of submm galaxies to constrain 
our best-fit model. At 850 iim, due to the K-correction, the 
source count is sensitive to the evolution of IR galaxies in 
a wide redshift range and out to very high redshift s. We 
showed that there is a degeneracy between the rate by which 
the characteristic luminosity of IR galaxies should increase 
to reproduce the source count and the maximum redshift 
out to which this increase should be continued; we resolve 
this degeneracy by requiring the model to reproduce the ob- 
served redshift distribution of submm galaxies. Moreover, 
we showed that our model requires a colour evolution to- 
wards cooler typical dust temperatures at higher redshifts. 
The employed colou r evolution is similar to that used by 
IValiante et al.l (120091 ) however, our best-fit model predicts a 
somewhat stronger colour evolution than that proposed by 
these authors. 

Another important feature of our model is that the best- 
fit is obtained using pure luminosity evolution but mildly 
evolving high-luminosity and low-luminosity slopes in the 
LF. Since high-luminosity sources are rare, the evolution 
in the high-luminosity slope is of little consequence. How- 
ever, the evolution of the low-lumuninosity slope affects large 
numbers of galaxies and if confirmed, this effect must have 
a physical origin, which can be addressed using numerical 
simulations of the evolution of the galaxy population, as well 
as with a combination of deep Herschel and optical imaging. 

The 850 /im-constrained best-fit model is consistent 
with observed 1100 /im and 70 /im source counts, which 
confirm the consistency of the implied colour-lumionosity- 
redshift distribution at both low and high redshifts. How- 
ever, this model under-produces the observed source counts 
at intermediate wavelengths, namely at 160 /im and SPIRE 
bands. To resolve this issue, we used the observed data at 
different wavelengths to find best-fit models which can repro- 
duce their observed source counts. While the best-fit models 
constrained at 70 ^m and 850 ^m are consistent with each 
other and also 1100 /xm, the implied evolutions for models 
capable of reproducing observed counts at other wavelengths 
are too diverse to be reconciled in a single model; specifically 
they need too strong colour evolutions which contradicts 
850 fim observations. While the inconsistenc y at 160 /im 
has been reported in earlier works (Le Borg ne et al.l 12009 : 
IValiante et al.ll2009l ). we are the first to report it for 250 /im, 
350 jum and 500 fim. We showed that the source counts at 
these wavelengths can be reproduced consistently, by adopt- 
ing the best-fit model which produces correct 70 /im, 850 fj,m. 
and 1100/im source counts, together with a modification 
in SED templates which is equivalent to the existence of 
a cold population of dusty galaxies at low to intermedi- 
ate redshifts which are under-represented in IRAS data. 
Besides the fact that there is some^ observational evidence 
for the existenc e of su c h galaxies (iStickel et al. 1998 . 2000 : 



Chapman et~ai] |2002|; [Patris et al.l |2003|; iDennefeld et al I 



|Sti 

M 



20051 : ISaiina et al.|[200^ : I Amblard et al.ll2O10l 'l. this assump- 



tion is similar to what other models ha d to assume in order 
to reproduce adequate 160 /im sources (|Lagache et al.ll2003l . 
12004 IValiante et al.|[2009l '). 

It is important to keep in mind that phenomenologi- 
cal models like what we described in this paper, are mainly 
simple mathematical forms which relate different observa- 
tions consistently rather than being physical models with 
explanatory power. However, their performance at differ- 
ent wavelengths and the distribution of sources they require 
for different redshifts can be used as their main predictions 
which also could be used to test their validity. While we used 
the redshift distribution of submm galaxies to constrain our 
model, we note that the observed redshift distribution of 
other wavelengths, if avail able, are in agreem e nt with our 
best- fit model predictions ([Jacobs et al.l I2OIII : iBerta et al.l 

Additional information including tabulated data for dif- 
ferential and cumulative source counts at different wave- 
lengths and their redshift distributions is available at 
http : //www. strw.leidenuniv.nl/genesis/ 
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possible colour, luminosity and redshift ranges into very 
small bins, assuming that in each bin the related variable is 
not changing significantly. 

The finite number of SED models we are using au- 
tomatically splits the colour range into 64 bins between 
0.29 < 7?(60, 100) < 1.64 (see Section [23J. We also use 
logarithmically spaced bins to split the possible luminosity 
range of IO^Lq < L < 10^'^Lq into 100 bins in our calcu- 
lation. This logarithmic scale which makes the integration 
roughly insensitive to the number of luminosity bins, is 
chosen to cope with the exponential nature of luminosity 
function where faint objects are much more numerous 
than luminous ones. It is also worth mentioning that the 
source count calculation is not sensitive to the minimum 
or maximum luminosity which is used in integration, if the 
used luminosity range covers the important 10^" — 10^'^Lq 
luminosity range; for instance, using Lmin = 10^ Lq instead 
of Lmin = lO^Z/Q as the minimum possible luminosity, does 
not change any of the source count calculations we are 
presenting in this paper. 



As discussed in Section [2.41 the uniform distribution of 
galaxies in redshift space is assured by the algorithm we are 
using, independent of the size of redshift bins. But, for the 
precise calculation of the K-correction and evolution func- 
tions, we split the redshift range of < 2 < 8 using bin 
sizes equal to Az = 0.01. However we noted it is possible 
to use even bigger redshift bins (e.g. Az = 0.1) without any 
significant change in the results. 



APPENDIX A: SOME NUMERICAL DETAILS 

For each specific evolution model, the source count at a 
given flux threshold and wavelength can be calculated 
based on equation (|13p . where the integration should be 
performed over all possible luminosities, colours and red- 
shifts. As mentioned in Section [2. 41 we do this by splitting 



© 0000 RAS, MNRAS 000, 000-000 



